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Abstract 

Tracking calculations of charged particles in electromagnetic fields require in principle the simul- 
taneous solution of the equation of motion and of Maxwell's equations. In many tracking codes a 
simpler and more efficient approach is used: external fields like that of the accelerating structures 
are provided as field maps, generated in separate computations and for the calculation of self fields 
the model of a particle bunch in uniform motion is used. We describe how an externally computed 
wake function can be approximated by a table of Taylor coefficients and how the wake field kick 

1 can be calculated for the particle distribution in a tracking calculation. The integrated kick, rep- 

P-i 

resenting the effect of a distributed structure, is applied at a discrete time. As an example, we 
use our approach to calculate the emittance growth of a bunch in an undulator beam pipe due to 
resitive wall wake field effects. 
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INTRODUCTION 



Self consistent particle tracking needs the simultaneous solution of the equation of motion 
for a multi particle system (N ~ 10 3 ■ ■ • 10 9 ) and of Maxwell's equations. The effort is 
usually rather high, as the electro-magnetic field calcualation in time domain has to resolve 
the dimension of the bunch and of the surrounding geometry: a large volume has to be 
discretized with high resolution, and a good spatial resolution is related to a fine time step 

Q. 

Supposed the particle energy is high enough for ultra-relativistic approximations, the 
motion is approximately uniform and the betatron wavelength is large compared to ge- 
ometric dimensions as cavities, magnets and wake generating discontinuities. Then the 
electromagnetic problem can be split into independent sub-problems: external fields, wake 
field interactions and simplified self-interactions. The separation of external fields is possi- 
ble without these conditions, but the wake field approach is based on an idealised setup as 
sketched in Fig. [1] a discontinuity, either of geometry or of material properties or of both, 
is enclosed by semi-infinite beam pipes. Several generations of wake field codes have been 
developed to solve this problem {2-6]. The model of simplified self interactions calculates so 
called "space charge forces" for a particle distribution in uniform motion, using a Lorentz 
transformation and solving an electrostatic field js iO]. 

This report is about the application of wake fields in a tracking program. The effect of a 
distributed wakefield generating structure is replaced by a discrete kick at a discrete time, 
calculated from a wake function that is determined with help of a wake field code. The 
wake function has to be described by a flexible format. We propose a hybrid formulation: 
the dependency on the longitudinal coordinate s (between source and observer particle) is 
tabulated and the transverse dependency (offset of source and observer particle) is Taylor 
expanded. Therefore the wake function is stored in a table of Taylor coefficents. The point 
to point interaction between all particles is calculated by a convolution method. Therefore 
the line charge density and generalized line charge densities (the first and second transverse 
moments) are calculated on a grid by a binning and smoothing technique. 

Our method is implemented in the space charge code ASTRA Jl]]. As an example we 
calculate the emittance growth of a bunch in an undulator beam pipe due to resitive wall 
wake field effects. In this example the wake per length in a 30 m beam pipe is modeled 
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by about 60 diskrete kicks. The example considers a real focussing lattice, collective and 
individual transverse motion and simultaneously longitudinal and transverse effects (which 
cannot be separeted from each other). 



COORDINATE SYSTEM AND WAKE FUNCTION 



The particles are tracked in global Cartesian coordinates with (x,y, z,p x ,p y ,p z ), the 
components of location and momentum, while wake fields are calculation in a local coodinate 
system with (u,v,w,p u ,p v ,p w ) components. The coordinates of both systems are related 
by: 
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with n the particle index, f the vector to the origin of the wake geometry, e w the unity vector 
into nominal direction of particle motion and e u , e v the orthonormal transverse vectors. The 
discrete wake induced kicks are added to the particle momenta when the center of mass of 
the distribution passes the plane f Q + ue u + ve w . The setup is skeched in Fig. [TJ 

For simplicity we set r D = for the following. A source particle with charge q s on the 
trajectory r s (t) = u s e u + v s e v + cte w is followed by a observer particle with charge q Q and 
trajectory f (t) = u Q e u + v Q e v + (ct — s)e w . The source particle causes the electromagentic 
field E^ s \ B^ s \ The total change of momentum of the observer particle is 



Ap = q (E^ s \f (t),t) + ce w x B^ s \f (t),t)) dt. 

J — oo 

For the definition of the wake function we follow j?! with 



(3) 



Wf(u s ,V s ,U ,V , s) 



q q s 



-Ap. 



(4) 



The wake function is causal in the s coordinate, the transverse and longitudinal components 
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are related by the Panofsky-Wenzel-Theorem 
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and the longitudinal component is a harmonic function of transverse coordinates 



(5) 
(6) 
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(wf ■ e w ) = 0. 
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For the following we use scalar component functions 



' h u (u s ,v s ,u ,v ,s) ^ 



Wf[U s ,V s ,U ,V ,S) 
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h v (u s ,v s ,u ,v , s) 
\yh w (u s: v s: u ,v ,s) J 

The implementation of wake fields in Astra is based on a second order Taylor expansion of 
the longitudinal wake in the transverse coordinates 



hyj(U s , V s , Uq, Vq, 
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This approach fulfilles Eq. (jJJ). The transverse components are uniquely related to the 
longitudinale wake by causality and Panofsky-Wenzel-Theorem: 

h u (u s , v s , u , v , s) = h$ (s) + 2hll ( s ) u s + 2/123 ( s ) v * + 2h 33 ( s ) u o + 2/134 (s)v , (10) 
h v (u s ,v s ,u ,v ,s) = h^l(s) + 2h\l{s)u s + 2h 2 l \{s)v s + 2h^l(s)u - 2h^(s)v Q , (11) 



with the integrated coefficient functions 



h%{s) 



h a p(x)dx 



(12) 



Special cases for geometries with symmetry of revolution are the monopole and dipole wake. 
The monopole wake 

Wf(u 8 , v s , u , v , s) = -h o(s)e w (13) 
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is purely longitudinal and independent on offset parameters. The transverse part of the 
dipole wake depends linear on the offset of the source particle. Due to symmetry the coeffi- 
cient functions hi 3 (s) and h 2 i(s) are identical. Therefore the dipole wake functions is 

w f (u s , v s , u Q , v Q , s) = -(u s e u + v s e v )2h^(s) - (u s u a + v s v )e w 2h 13 (s) . (14) 



DISTRIBUTED SOURCE 

For an arbitrary distribution of source particles with charges q s ^ n , total charge qt, = ^ q s ,n 
and trajectories r s n (t) = u s ^ n e u + v s ^ n e v + (ct + w s>n )e w , the change of momentum of the 
observer particle is 



Ap= — Qs,nWf(u Stn , V s>n , U , V , S - W Stn ) . (15) 
n 

The wake potential 

W(u Q , v , s) = — Y] q s , n Wf(u Stn , v s nt u , v , s - w Sin ). (16) 
qb 

n 

characterises the shape dependent wake kick. It can be computed with electromagnetic 

u n n 

field calculation programs as [2|, 13|, |6J for continuous gaussian source distributions. Usu- 
ally the wake function is extrapolated from wake potential computations for small source 
distributions. 

The particle summation in Eq. ffT5"]) is replaced by a volume integration 



Ap = —J g(u,v,w)wf(u,v,u ,v , s — w)dudvdw (17) 

with g(u, v, w) the charge density. Formally this can be expressed with a Taylor expansion 
of the wake function as 

Ap = — J g(u, v, w) (wij mn (s — dudvdw . (18) 

ijmn 

Changing the order of summation and integration Ap can be written as sum of one dimen- 
sional convolution integrals 

A P = ~ ^ u o v o J hj{w)w ijm n{s - w)dw , (19) 



ijmn 
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with scalar functions 

\ij(w) = J g(u,v,w)u l v^dudv , (20) 

and Wij mn (s) described in Tab. HI Xqq(w) is the line charge density. In the following we call 
Kj(w) generalized line charge densities. 

NUMERICAL REALIZATION 

For the numerical kick calculations three problems have to be solved: the calculation 
of generalized one dimensional density functions \j(s), the representation of coefficient 
functions h a p{s) and the convolution of these functions. 

Binning and Smoothing 

A continuization based on a binning and smoothing technique is used to convert the 
particle set with longitudinal positions {w n } and discrete generalized charges {q Stn u l n vl} 
to continuous generalized line charge densities Ajj(s). Conventional binning distributes the 
particle charges to N^ n equi-spaced bins ranging over the whole bunch length. The disad- 
vantages are a spacial resolution independent on particle density, and the bin boundaries 
are determined by extreme particles (first and last particle). To achive more flexibility, we 
use bins that fulfill requirements for the length per bin as well as for the charge per bin, and 
we utilize a multiple binning with N su b different conditions for the first and last bin. 

Therefore the bunch is split into Nbi n N su b sub-bins. All inner bins are composed by N su b 
consecutive sub-bins, the outer bins may combine less sub-bins, as sketched in Fig. [2J By 
doing this, N su b different binnings are generated. 

The sub-binning is controlled by the linear combination f(w) = W eq f w (w) + (1—W eq )f q (w) 
of two auxiliary functions f w (w) and f q {w). These functions increase along the bunch 
monotonously from zero to one. Function f w (w) is proportional to length while f q (w) 
increases with the integrated bunch charge. The boundaries bk of the sub-bins are the 
solutions of fipk) = k / (N b i n N su b) . Setting W eq to one causes equi-spaced sub-bins, zero 
leads to equi-charged sub-bins and a middle value generates a mixture of both conditions. 
The same bins and sub-bins are used for all generalized charges. 

Result of the binning procedure are Nbi n N su b + N su b — 1 different bins, each characterized 
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by the bin center b Cyk = {b k + b k+Nsub )/2, the bin width b w , k = (b k+Naub - b k ) and the bin 
weights (the sums of generalized charges) 



Nsub 

nebm fe 

The smoothing procedure replaces all bins by continuous functions with certain centers, 
widths and weights and samples the sum of these functions on an equidistant grid: 

Kj,n = Xij(nAw) , (21) 
A*» = Y^T^SC^) . (22) 

. "w,k Vw,k 

The ASTRA implementation uses rectangular, triangular and gaussian smoothing functions: 

\ 1 |z| < 0.5 
S r (x) = < 

I otherwise 

S t (x) = max(0, 1 — |x|) , 
with p ~ 1, a control parameter. The grid density is determined by the smallest bin width. 



Representation of Coefficient Functions 

Each coefficient function h a p(s) is described by coefficients R a /3, L a p,C aj 3 and limited 
auxiliary functions q a p{s), p a p(s). For the following we skip the index. The representation 
is 

Hs) = q(s) + ^ + RcS(s) - [Lc5(s) + p(s)} , (23) 

with S(s) the dirac function, $(s) the step function and C = C if C ^ 0. The term with 1/C 
is skipped for C — 0. The representation is not unique, but very flexible: every coeffcient 
R, L, C might vanish and the auxiliary functions migth be identical to zero. The auxiliary 
functions are polygons as sketched in Fig. [3j The table description of auxiliary function 
h a p{s) is listed in Tab. Ull The entry "10a + /3" describes the subscript of the auxiliary 
function. For a complete wake field representation all sub-tables of non vanishing coefficient 
functions are stacked to one table 

Table = stack 
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N t ,Tablel,Table2. . .j 



with Nt the number of sub-tables. The order of sub-tables is arbitrary, vanishing coefficients 
need no representation. 



Convolution 

According to Eqs. ( !T9|20|) we need convolution integrals of generalized charges with co- 
efficient functions and integrated coefficient functions: 



H(s) = J x{x)h{s ~ x)dx • 

H®(s) = [ X(x)h^(s-x)dx 



For simplicity we skipped all lower indices. Using the representation Eq. ( 123]) and auxiliary 
functions 

A(s) = q(s) + <S>(s)/C , 
B(s) = -cp(s) , 
Q(s) = / q(s)ds , 

J — oo 

C(s) = -Q(s) + cp(s) - (s/C + Rc)$(s) , 

the convolutions can be rewitten as 

H(s) = RcX(x) - Lc 2 \'(x) + J \{x)A{s - x)dx + J \'{x)B{s - x)dx , (24) 

H®(s) = Lc 2 X(x) + J \{x)C(s-x)dx . (25) 

The integrals are solved piecewise for the generalized charged densities Eq. (1221) and polyg- 
onal coeffcient functions. 

EXAMPLE 

The example of particle motion in an undulator vacuum chamber with resisitve wall wake 



fields is based on parameters used in 11], see Tab. IIII1 For the calculation with ASTRA, a 
round stainless steel beam pipe with radius r& = 4.75 mm is used to cause significant trans- 
verse effects. The undulator chamber of the real facility (FLASH) is made from aluminium. 
For the numerical simulation 63 quadrupoles of the FODO lattice are considered and the 
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resistive wall wake of the 30 m beam pipe is simulated by 62 discrete wake kicks (each for 
AL = Lfodo/2) applied in the middle of the drifts between the quadrupoles. 

The longitudinal wake per length of a round resistive beam pipe is approximated by 

Wf(u s ,v s ,u ,v , s) -e w = w r {s) ^1 + 2 j , (26) 

with 

1 r 

w r (s) = / Z r {uS) exp(jus/c)duj , 

27T J 

Z r (u>) = ^(1+^(^/2)"' , 
2irr b 

Z s (u) = y/jufJ,/K(u) , 

\i the permeability, e the permittivity of vaccum and k(uj) the frequency dependent conduc- 
tivity of the pipe [l^]. therefore only three coefficient functions 

^oo (s) = -ALw r (s) , 
ALw r (s) 



hi 3 (s) 



, , . ALwr(s) 
h 2 As) = - 



have to be calcualted that are all proportional to w r (s), the longitudinal resistive wall wake 
per length. In the wake table the coefficients N p , R, L und C are set to zero (for a/3 = 00, 
13 and 24), only the auxiliary functions Woo(s), Uis(s) and ^24(3) are used. In Figs. H] 
the emittance growth due to dipole wake forces is shown along the undulator and compared 



with the analytic estimation from [U| . The numerical result agrees well with the estimation, 
but the dipole model is incomplete. Fig. \5\ shows the emittance growth due to monopole 
and dipole wake fields. The longitudinal wake causes an energy spread of abaut 2.4 MeV 
after 30 meters. For off axis particles with that energy spread, the quadrupole focussing is 
significantly altered and the transverse emittance is further increased. 
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TABLE I: Wij mn (s) function. 
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N q 




R /(Us) 


L /(Us 2 ) 


C /(1/U) 


10a + /3 


s q ,\ /m 


q(s q ,i) /u 


Sq,2 /m 


9(s g , 2 ) /U 






Sq,N q /m 


^ (•Sg,AT„) /U 


s p ,i /m 


p(a Pjl ) /(Us) 


s P ,2 /m 


P(«P,2) /(Us) 




... /(Us) 


s P ,N q /m 


P(S P ,7V P ) /(Us) 



TABLE II: Sub-table of auxiliary function h a p(s). The unit U is V/As for a/3 = 00, V/(Asm) for 
a/3 = 01 • • • 04 and V/(Asm 2 ) for all other coefficients. 



variable 


units 


values 


beam energy 


£o 


GeV 


1 


bunch charge 


Qo 


nC 


1 


rms bunch length 




/im 


50.0 


normalized emittance in the undulator 


e n 


/im 


2.0 


undulator length 


L 


m 


30.0 


undulator gap 


9 


mm 


12.0 


pipe thickness (minimum) 


t 


mm 


1.25 


beam optics beta function 




m 


3 


FODO period length 


Lfodo 


m 


0.96 


Length of FODO quad 


Lfodo 


mm 


136.5 



TABLE III: Parameters for wake field calculation in undulator as asumed for TTF2. 
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FIG. 1: Source particle (red) and observer particle (green) in a cavity between semi infinite beam 
pipes. 
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FIG. 2: Binning and sub-binning for Nun = 4 and N su b = 3. 
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FIG. 3: Auxiliary function q(s) of poligonial shape, with si < S2 < ■■■ < %. 
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FIG. 4: Emittance growth along the undulator for dipole wakes in a stainless steel vacuum chamber 
with a radius of 4.75 mm: numerical simulation (solid line); analytical solution (dashed line). 
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FIG. 5: Emittance growth along the undulator for monopole and dipole wakes in a stainless steel 
vacuum chamber with a radius of 4.75 mm: numerical simulation (solid line); analytical solution 
(dashed line). 
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